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[i] A standing shock wave behind the Moon was predicted by Michel (1967) but never 
observed nor simulated. We use ID hybrid code in order to simulate the collapse of the 
plasma-free cavity behind the Moon and for the first time to model the formation of this 
shock. Starting immediately downstream of the obstacle we consider the evolution of 
plasma expansion into the cavity in the frame of reference moving along with the solar 
wind. Well-known effects as electric charging of the cavity affecting the plasma flow and 
counterstreaming ion beams in the wake are reproduced. Near the apex of the inner Mach 
cone where the plasma flows from the opposite sides of the obstacle meet, a shock wave 
arises. We expect the shock to be produced at periods of high electron temperature solar 
wind streams (T ; c T e ~ 100 eV). The shock is produced by the interaction of oppositely 
directed proton beams in the plane containing solar wind velocity and interplanetary 
magnetic field vectors. In the direction across the magnetic field and the solar wind 
velocity, the shock results from the interaction of the plasma flow with the region of the 
enhanced magnetic field inside the cavity that plays the role of the magnetic barrier. The 
appearance of the standing shock wave is expected at the distance of ~7R M downstream of 
the Moon. 

Citation: Israelevich, P., and L. Ofman (2012), Hybrid simulation of the shock wave trailing the Moon, J. Geophys. Res., 117, 
A08223, doi: 10. 1029/201 1JA017358. 


1. Introduction 

[ 2 ] An idealized picture of the interaction between 
supersonic and super-Alfvenic solar wind flow and non- 
magnetized celestial bodies without an atmosphere such as 
the Moon, asteroids, etc. (see Figure 1) can be outlined as 
follows [e.g., Michel, 1968]. The solar wind particles are 
absorbed by the surface when they impact the celestial body 
whereas the velocity of the bypassing particles remains 
unchanged. As a result, a cavity free of solar wind plasma is 
formed behind the obstacle. Because of thermal dispersion 
of velocities of the impacting particles population, some of 
the particles behind the limb have a velocity component 
directed toward the axis of symmetry of the cavity. These 
particles penetrate the geometric shadow of the celestial 
body. Since the bulk velocity of the solar wind is larger than 
the thermal velocity, the depth of penetration is proportional 
to the distance from the limb and the cavity free of plasma 
takes a conical shape. Considering the frame of reference 
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moving with the solar wind velocity, the process of refilling 
of the region of geometrical shadow with plasma is equiva- 
lent to the plasma expansion in vacuum [Samir et al., 1983]. 
Since the thermal velocity of electrons is much higher than 
the solar wind speed therefore, at the very initial stage fast 
electrons enter the whole region of the geometrical shadow 
charging it negatively. As soon as the negative potential 
reaches ~ TJe (where T e is the electron temperature) further 
charge separation stops, and electrons continue to fill the 
cavity along with solar protons that are accelerated by the 
negative potential in the cavity to the ion-acoustic velocity 

c s ~ y/T The cavity is restricted by the inner Mach cone of 

the angle tan (3 = cjv sw , where v sw is the solar wind speed 
[Michel, 1968]. Plasma expansion into the cavity is accom- 
panied by rarefaction wave propagation away from the cavity 
axis. Since the solar wind is supersonic the outer Mach cone 
restricts the region occupied by the rarefaction wave and the 
solar wind plasma flow outside this cone remains 
unperturbed. 

[ 3 ] The presence of the interplanetary magnetic field 
(IMF) introduces asymmetry in this description. If we con- 
sider the coordinate system with X-axis along the solar wind 
outflow velocity and 7-axis defined in a way that the IMF 
vector lays in the XT-plane, the plasma flow along the 
magnetic field fills the cavity. If the R r -component of the 
magnetic field is present and the magnetic field lines are 
inclined to the 7-axis at an angle a = arctan (B x /B y ), then the 
bulk plasma velocity directed toward the wake axis is dif- 
ferent at the opposite sides of the wake as shown by vector 
diagrams in Figure 1. The inner Mach cone is not 
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Figure 1 . A diagram showing the formation of the lunar wake. The arrows show the directions of the 
solar wind velocity and the interplanetary magnetic field (IMF). Thick solid lines show the positions of 
the boundaries, blue, yellow and red colors correspond to the cavity free of plasma, rarefaction region 
and the region of enhanced plasma density, respectively. Vector diagrams explain the plasma entrance into 
the cavity from two opposite sides. Bottom panel shows the wake structure for the case when the IMF is 
aligned along the solar wind velocity. 


symmetric: the angles at the apex of the cone are not equal 
(3\ ± /? 2 - Simple geometric consideration shows that 
c coso; 

tan P 12 = 7 , and the length of the cavity is L = 

' v sw ± c s sin a 

M 2 — sin 2 a 

2 R m — , where R M is the Moon radius and M s = 

M s cos a 

v sw lc s is the Mach number. For the case B x = 0, L = 2 R m M s . 
The cavity length doubles for a = 60° and grows rapidly 
with further increase of a. For case of IMF parallel to the 
solar wind velocity (a = 90°), the cavity becomes an infinite 
cylinder as shown in the bottom panel of Figure 1 . 

[ 4 ] Across the magnetic field in the Z direction, plasma 
moves along with the frozen-in magnetic field lines. This 
flow results in magnetic field compression inside the cavity 
and magnetic field decrease in the rarefaction region 
between the inner and outer Mach cones. 

[5] These main features of the solar wind interaction with 
the Moon were observed in the very first measurements 


by Explorer 35 satellite [Colburn et ai, 1967; Lyon et al., 
1967; Ness et al, 1967] and later confirmed by the Wind 
[Ogilvie et al., 1996; Owen et al., 1996; Bosqued et al., 
1996], LRO [Halekas et al, 2005], Kaguya [Nishino et al., 
2009], and Artemis [Halekas et al., 2011] spacecraft. They 
were reproduced in laboratory experiments as well 
[Kristoferson, 1969; Podgorny and Andrijanov, 1978; 
Podgorny et al, 1980], and in numerical simulations 
[Lipatov 1976, 2002; Farrell et al., 1998; Birch and 
Chapman, 2001; Kallio 2005; Kimura and Nakagawa, 
2008; Lipatov et al, 2011; Wang et al., 2011, Wiehle et al., 
2011; Holmstrom et al., 2012]. 

[6] However, there is another feature expected for the 
solar wind interaction with the Moon that was not observed 
or modeled yet. Michel [1967, 1968] pointed out that a 
consequence of solar plasma flow into the cavity would be 
the appearance of a shock wave at the place where the col- 
lapse of the cavity is halted. The orbits of previous near- 
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Moon spacecrafts were too low in order to reach the region 
of expected shock and it was not observed experimentally. 

[ 7 ] One of the possible explanations why the shock wave 
was not observed in past numerical simulations is as follows. 
The typical electron temperature used in these simulations 
was T e ~ T, ~ 10 eV [. Lipatov 1976, 2002; Farrell et ah, 
1998; Birch and Chapman, 2001; Kallio 2005; Kimura 
and Nakagawa, 2008; Lipatov et al, 2011; Wang et al., 
2011, Wiehle et al., 2011; Holmstrom et al., 2012], In the 
present study, we model for the first time the formation of 
this shock for the case T e » T, in hot solar wind streams 
(T e ~ 100 eV) using ID hybrid code. 

[8] Near the apex of the inner Mach cone two opposing 
flows of plasma meet and their relative velocity is approxi- 
mately 2 c s . Therefore it is natural to expect that a standing 
shock will arise at this location as shown in Figure 1. In case 
of the magnetic field parallel to the solar wind flow (the 
cavity is not collapsed as shown in the bottom panel of 
Figure 1) the standing shock wave may still exist. If the 
transition from collapsing cavity to a region of enhanced 
magnetic field with constant cross-section occurs as dis- 
continuous change in local physical quantities then the 
appearance for the shock wave is evident. A smoothing out 
of the transition between the two limiting shapes does not 
eliminate the shock wave except in the immediate vicinity of 
the boundary where it becomes a weak discontinuity that 
evolves in a shock wave [Michel, 1968; Wolf, 1968]. In this 
paper, we perform for the first time hybrid numerical simu- 
lation of the shock wave formation behind the Moon. 


periods, while neglecting the small temporal- and spatial- 
scales of the electron kinetic motions. In the hybrid model, 
the three components of order million particle velocities 
were used to calculate the currents, and the fields on the ID 
grid. In our study the grid consisted of 512 cells, and the 
total length of our simulated region was 139 c/w p , where w p 
is the plasma frequency, and c is the speed of light. The 
required number of particles per cell was determined by the 
required limitation on the overall statistical noise, and we 
have used up to 1600 particles per cell to reduce statistical 
noise to acceptable levels (i.e., the random noise in velocity 
space is much less than the thermal distribution). Note, that 
each numerical particle represents large number of real par- 
ticles, with the conversion determined by the density 
normalization. 

[ 12 ] The following equations of motion are solved for the 
ions: 



(1) 


, . d\, ( Vi x B 

M — — = e E H 

at V c 


(2) 


where M is the proton mass and e is the electron charge. The 
electron momentum equation is solved by neglecting the 
electron inertia 


—n e m e y e = 0 = -en e 
at 


v ( xB\ „ 

E H ) - X7p e 


(3) 


2. Model 

[ 9 ] We follow the simple method used by Farrell et al. 
[1998] and Birch and Chapman [2001]. Let us consider 
the frame of reference moving along the X-axis with the 
solar wind velocity v sw . In this frame, the solar wind plasma 
expands into the plasma-free region with initial shape 
corresponding to the cross-section of the obstacle. Solution 
of the time-dependent 2D fluid problem is f{t, v, z). Making 
the substitution x = v sw t, one returns to the obstacle frame of 
reference and obtains an approximate solution of the 3D 
problem of steady state flow around an obstacle /(x, y, z). 
Similarly, ID problem of plasma expansion in vacuum 
considered in this study yields the approximate solution 
fit, y) — > f{x, v) for 2D flow around an infinite plate of 
width 2y 0 with an axis parallel to Z. 

[ 10 ] We consider supersonic and super Alfvenic plasma 
flow around an infinite cylinder. The coordinate system is as 
follows: X-axis is along the solar wind velocity, Z-axis is 
parallel to the cylinder axis (d/dz = 0), and T-axis completes 
the system. Solutions are obtained for two cases (a) expan- 
sion along the magnetic field (B 0 = (0, B 0 0)), and (b) 
expansion across the magnetic field (B = (0, 0, B 0 )). 

[ 11 ] We use the ID hybrid code, based on the code 
developed by Winske and Omidi [1993], and later modified 
by A. Vinas (private communication, 1995) and Ofman et al. 
[2001, 2009]. In the hybrid model the ions are represented as 
particles, neglecting collisions, while the electrons are 
described as a finite temperature massless fluid that main- 
tains the quasineutrality of the plasma. This numerical 
modeling method allows one to resolve the ion dynamics 
and to integrate the equations over many ion cyclotron 


where the electron pressure p e = n e T e is used for closure, and 
quasi-neutrality implies n = n e = n„ where n k is the number 
density of electrons, and protons, respectively. Thus, the 
collisionless Ohm's law is obtained. We use the isothermal 
energy equation for electrons (7 = 1, T= T 0 = Const.), hence 
p e is proportional to the number density n. 

[ 13 ] The above equations are supplemented with Max- 

471 1 SB 

well's equations V x B = — J and V x E = — . The 

c c at 

field solutions are obtained on the ID grid, and the particle 
(ion) equations of motions are solved as they respond to the 
fields at each time step. Three components of magnetic field 
and velocity vectors are included in the model. The numeri- 
cal method has been tested and used successfully in many 
studies. 

[ 14 ] The particle and field equations are integrated in time 
using the Rational Runge-Kutta (RRK) method [Wambecq, 
1978], whereas the spatial derivatives are calculated by pseu- 
dospectral FFT method. The hybrid model allows computing 
the self-consistent evolution of the velocity distribution of the 
ions that includes the nonlinear effects of wave-particle inter- 
actions without additional assumptions. Moreover, the hybrid 
model is well suited to describe the nonlinear saturated state of 
the plasma. 

[ 15 ] Typical parameters of the solar wind plasma at 1 AU 
are chosen as the initial condition: «o = 5 cnf 3 , Bo = 10 nT, 
Tj = 5 eV and the size of the obstacle 2y 0 = 3500 km is the 
Moon’s diameter. 

[ 16 ] Since the necessary condition of ionacoustic instabil- 
ity development is T, c T e , we chose the electron temperature 
to be much higher than ion temperature: T e = 100 eV. The 
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Figure 2. Distribution of the (a) number density, (b) electric potential, and (c) plasma velocity compo- 
nent V y in the wake. Solid lines denote the plasma expansion with ion-acoustic velocity c s . Expansion with 
V = 220 km/s is shown in Figure 2c. 
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Figure 3. (a) Plasma density and (b) electric potential distribution inside the cavity. Solid lines show 
expansion with velocities V = 150 km/s (cavity boundary) and V = 220 km/s. 


electron temperature in the solar wind can reach occasionally 
40 eV [Newbury et al., 1998]. Moreover, the solar wind 
electrons distribution function is not Maxwellian [see, e.g., 
Feldman et al., 1975]. It rather consists from cold core 


electrons with the temperature 7], ~ T, ~ 1 0 eV, and hotter 
halo (~7%) with the temperature ~50-100 eV. Sittler and 
Burlaga [1998] found that within solar wind magnetic 
clouds T e » Tj, and T e / 7) ~ 7, on average. Thus, our value of 
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Figure 4. Distribution of protons in phase space (v — v v ) at 9 distances along the wake. 


T e represents the extreme conditions that guarantee the for- 
mation of the shock. However, we found that the shock can 
form at lower T e as well (see below). 


[n] The obstacle is centered aty = 0 at the time t=0 (i.e., at 
the distance x = 0). The plasma density is n 0 for [y| > y 0 , and 
0.1/?o behind the obstacle for [y| <_v 0 (zero density cannot be 
used for numerical reasons). All other parameters (velocity, 



Figure 5. Proton velocity distribution functions at X— 16,000 km for Y= —300 km (blue), Y = 0 (green), 
and Y = 300 km (red line). 
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Figure 6. Dependencies of proton velocity distribution 
function on the distance in the wake X for the wake sections 
shown by dashed lines in Figure 2a. From top to bottom: 
7 = 2800 km, Y= 1000 km, Y= 300 km, and 7=0. 

temperature, magnetic field) are constant along the 7-axis at 
t = 0 (x = 0). The solar wind velocity v sw = 500 km/s is not used 
explicitly in the simulations but serves for conversion of the 
time coordinate to the distance along the X-axis: x = v sw t. 

3. Results 

3.1. Expansion Along the Magnetic Field (B x = B z = 0) 

[is] Figure 2a shows the color-coded plasma density 
behind the obstacle. Simulation reproduces both the col- 
lapsing cavity and expanding rarefaction region restricted by 
the outer Mach cone. The cavity collapses at X = 6000 km. 
Starting from this distance, accumulation of plasma occurs 
near the center of the wake. The region of the enhanced 
plasma density expands with the distance reaching 
^ 1 500 km atX= 16,000 km. Rarefaction continues between 
the outer boundary of this region and outer Mach cone. 
Narrow (~300 km) minimum of the density is observed in 
the closest vicinity of the wake axis (7= 0) starting from the 
distance of X ~ 10,000 km. 

[ 19 ] For the neutral gas expansion in vacuum, the apex 
angles of the inner and outer Mach cones are the same and 

equal to tan f3= VJV (where V s ~ and V are sound and 

flow velocities, respectively). Simulation shows that this is 
not the case for plasma expansion in vacuum. Thick solid 
lines in Figure 2a denote the inner and outer Mach cones for 
ion-acoustic velocity tan (3 = c Jv sw (c x = 98 km/s). Figure 2a 


shows that the cavity collapse occurs with the velocity 
145 km/s, which is higher than c s . Rarefaction wave forming 
the outer Mach cone propagates with 88 km/s, i.e., slower 
than c s . The asymmetry of perturbation propagation toward 
and away from the cavity is associated with the appearance 
of the electric field in the rarefaction region. 

[ 20 ] Figure 2b shows the distribution of the electric 
potential in the wake of the obstacle. The potential is 
strongly negative inside the cavity reaching —260 V. This 
negative potential results from the strong gradient of the 
number density at the boundary of the cavity, as it follows 
from the Ohm’s law (3): n = n 0 exp(ecp/T e ) or cp = ln/;/«o • 

[ 21 ] Figure 2c shows the distribution of the velocity 
component v v . Solid lines correspond to the expansion with 
ion-acoustic velocity c s and with the velocity V — 220 km/s. 
Thus, behind the solid body a conical region is formed (with 
the apex angle tan a = V/v sw ), where the number density 
remains the same as at the initial state. Outside this cone, 
protons are accelerated by the floating potential that appears 
inside the cavity. This results in the increase of the number 
density and the change of the electric potential. Further refill 
of the cavity occurs under the action of the electric potential 
at the boundary of the cavity. 

[ 22 ] Figure 3 shows the density (Figure 3a) and the elec- 
tric potential distribution (Figure 3b) inside the cavity. Due 
to the protons accelerated by the floating potential, the 
density inside the cavity (but outside the cone determined by 
the velocity V= 220 km/s) increases up to 0.8 cm' 3 , and the 
negative potential at the cavity axis drops from —260 V to 
— 180 V. Electric potential of the cavity boundary is 
U b = —120 V, so the velocity of the cavity collapse should 


be V = =150 km/s which is in good agreement with 

the simulation results. 

[ 23 ] Evolution of the proton distribution in phase space 
(y,v v ) is shown in Figure 4. Each panel represents the 
dependence of the distribution function on 7 (horizontal 
axes). Vertical axes are particle velocity component v y and 
the number of particles within the velocity range (v,,, v y + dv v ) 
is color-coded (the distribution functions are normalized in 


such a way that n = J f(v)dv). The results fort 


= 1 s, 5 s, 9 s, 


13 s, 17 s, 21 s, 25 s, 29 s and 32 s (which correspond to the 
distances X along the flow direction of 500 km, 2500 km, 
4500 km, 6500 km, 8500 km, 10,500 km, 12,500 km, 14,500 
km and 16,000 km, respectively). The average velocity starts 
to increase at the outer Mach cone, and grows linearly till the 
boundary of the inner Mach cone. Inside the cavity the 
populations of protons accelerated by the floating potential to 
240 km/s can be seen (X= 2500 andX= 4500 km). Starting at 
X = 6500 km two counterstreaming proton beams are 
observed near the axis of the wake (7=0). Also, the beams of 
high speed protons coming from the opposite sides of the 
obstacle can be seen in the rarefaction region starting from 
X = 8500 km. The proton beams in the rarefaction region 
practically do not interact with the ambient plasma streaming 
toward the wake axis and can be seen until the end of the 
simulation (X = 16,000 km). On the contrary, oppositely 
directed proton beams strongly interact in the vicinity of the 
axis. AtX= 8500 km, two beams are clearly seen, whereas at 
X = 10,500 km the beam structure disappears and at 
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Figure 7. Profiles of plasma density, macroscopic velocity component V y , and electric potential at 
X= 16,000 km. 
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X = 12,500 regions of standing plasma appear on both sides 
of the axis. These regions increase with the distance reaching 
±750 km from the axis. They correspond to the increase 
of the plasma density in the middle of the wake shown in 
Figure la. 

[ 24 ] In a very narrow region, Y < ±150 km, two beams 
of protons with velocities ^90 km/s continue to exist till 
the end of simulation. Here, narrow minima of the plasma 
density (Figure la) and electric potential (Figure lb) are 
observed. 

[ 25 ] Figure 5 shows the velocity distribution functions at 
X= 16,000 km at three different positions along 7: —300 km, 
0, and 300 km. Two beams (~±90 km/s) exist at Y= 0, and 
at the distance of 300 km two slower beams (±20 km/s 
and +40 km/s) almost overlap forming nearly Maxwellian 
velocity distribution. 

[ 26 ] Change of the velocity distribution function along 
Y-axis for given Y positions (Y= 2800 km, 1000 km, 300 km, 
and 0, denoted by dashed lines in Figure 1) is shown in 
Figure 6. For Y= 2800 km, the velocity distribution function 
remains unperturbed until X= 6000 km, when the outer Mach 
cone is crossed. Further down the wake, the velocity distri- 
bution function becomes shifted because of the plasma flow 
toward the wake axis. High-speed protons from the opposite 
sides of the wake are seen starting with X ~ 10,000 km. The 
line along Y = 1000 km starts inside the cavity and low- 
density plasma at rest is observed until X = 2000 km. Fast 
protons accelerated by the floating potential start to appear at 
this distance. The main injection of protons into the cavity 
starts at X = 3000 km. The flow velocity reduces from 
140 km/s to 105 km/s at X= 8000 km and further remains 
constant. The flux of protons from the opposite side of the 
wake appears at X= 5000 km and two proton beams coexist 
for X> 5000 km without any significant interaction. 


[ 27 ] For the line Y = 300 km, accelerated protons with 
v = 220 km/s appear at X= 3000 km, and the main plasma 
influx occurs atX= 5500 km. The oppositely directed beam 
of accelerated protons from the other side of the wake 
appears at X = 4000 km. The velocity of particles in this 
beam decreases with the distance, and two beams exhibiting 
some interaction exist until X= 1 1,000 km. At this distance, 
two-beam structure disappears, the plasma is at rest, and its 
density has increased. 

[ 28 ] Along the wake axis, 7=0, protons accelerated from 
both sides appear at X= 4000 km. At X= 6000 km the cavity 
collapses and the two beam distribution function becomes 
apparent. The relative velocity of the two beams reduced 
somewhat until X = 8000 km, and then remains practically 
unchanged exhibiting only quasiperiodic variations. 

[ 29 ] The bottom row in Figure 4 (X= 12,500 km, 14,500 km 
and 16,000 km) show the existence of a discontinuity sepa- 
rating the region of plasma flow toward the wake center 
and the region of the plasma at rest with enhanced number 
density. Vertical dashed lines denote the location of the 
discontinuity. The macroscopic plasma parameters also have 
discontinuity at that location, as shown in Figure 7. Here the 
profiles of the density, macroscopic velocity, and electric 
potential along the 7-axis are given for X= 16,000 km. The 
dashed lines show the discontinuities positions on both sides 
of the wake axis. 

[ 30 ] The simulated discontinuity has nonzero velocity 
component normal to the discontinuity surface and the 
plasma density and parallel velocity jump across the dis- 
continuity. The discontinuity propagates with the velocity 
41 km/s away from the tail axis. Behind the discontinuity the 
plasma is at rest. The incoming plasma velocity is 60 km/s in 
front of the discontinuity. Thus, in the discontinuity rest 
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Figure 8. Distribution of the (a) number density, (b) magnetic field, and (c) plasma velocity component V y in 
the wake for the case of expansion across the magnetic field. Solid lines denote the plasma expansion with ion- 
acoustic velocity c s . Expansion with fast magnetosonic velocity Vj- = 320 km/s is shown in Figures 8b and 8c. 
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Figure 9. Distribution of protons in phase space (y — v y ) at 9 distances along the wake for the case of 
expansion across the magnetic field. 


frame the plasma velocity is larger than c s on the external 
side, and less than c s behind the discontinuity. Therefore, 
this is collisionless slow shock wave [e.g., Swift, 1983]. 

3.2. Expansion Across the Magnetic Field (B x = B y = 0) 

[ 31 ] For the case of plasma expansion across the magnetic 
field, the cavity containing the magnetic field is compressed 
by the magnetized plasma flow. Because of magnetic flux 
conservation the magnetic field in plasma-free cavity 
increases. The compression ceases when the magnetic pres- 
sure in the cavity reaches the value of the plasma thermal 
pressure in the unperturbed plasma flow. This is illustrated 
in Figure 8a showing the density distribution behind the 
obstacle. Contrary to the expansion along the magnetic field, 
the compression of the cavity as well as rarefaction wave 
expansion occur with the same velocity approximately equal 
to the ion-acoustic speed c s , thus it can be described as the 
slow magnetosonic wave propagation. The pressure balance 
is reached aXX— 5000 km where the cavity shape is changed 
and for X > 5000 km, the cavity cross-section remains con- 
stant. The cavity is surrounded by the region of enhanced 
plasma density that starts at the transition from collapsing to 
a constant cross-section as was predicted by Michel [1968]. 
The region of enhanced plasma density expands with the 
speed c s . 

[ 32 ] Magnetic field distribution and the velocity compo- 
nent V y of plasma flow behind the obstacle are shown in 
Figures 8b and 8c. If the initial plasma density was zero 
immediately behind the obstacle, the magnetic field 


compression inside the cavity would propagate as an elec- 
tromagnetic wave in vacuum with the speed of light c. 
However, the plasma density behind the obstacle was chosen 
0.1 h 0 rather than zero (for stability of the numerical solu- 
tion), and the compression inside the cavity occurs with fast 
magnetosonic speed Vf = \J Kj + c\ = 320 km/s as it is seen 
in the figure. 

[ 33 ] Plasma flow toward the wake axis exists inside the 
rarefaction region. It vanishes rapidly at the boundary of the 
region of the enhanced plasma density. 

[ 34 ] Proton velocity distribution functions in phase space 
( y , Vy) are shown in Figure 9 for X = 500 km, 2500 km, 
4500 km, 6500 km, 8500 km, 10,500 km, 12,500 km, 
14,500 km and 16,000 km. 

[ 35 ] In contrast with the expansion along the magnetic 
field, there are no regions of two counterstreaming proton 
beams since the magnetic field prevents access of particles 
from one side of the obstacle to another side. As the mag- 
netic field pressure in the cavity becomes stronger the 
plasma flow toward the axis of the wake decelerates (as seen 
for X= 4500 km). For the non-collapsing part of the cavity, 
the velocity at the cavity boundary is zero, the density 
increases, and the region of the plasma at rest with enhanced 
density expands with the distance from the obstacle. Thus 
we observe the typical pattern of plasma flow against a 
magnetic barrier. The velocity exhibits a discontinuity at the 
boundary of the dense plasma region so that the outer 
boundary of this region forms a shock wave. 
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Figure 10. Profiles of plasma density, macroscopic veloc- 
ity component V v , and magnetic field at X = 16,000 km for 
the case of expansion across the magnetic field. 

[36] Figure 10 shows the profiles of the density, macro- 
scopic velocity and magnetic field along the 7-axis at 
X = 16,000 km. The quasi-perpendicular shock is observed 
in the case of plasma expansion across the magnetic field. 

4. Conclusion 

[ 37 ] The ID hybrid simulation shows for the first time the 
appearance of the shock wave trailing the Moon in the dis- 
tant wake, in agreement with past predictions [Michel, 
1967]. The shock is produced by the interaction of oppo- 
sitely directed proton beams in the plane containing solar 
wind velocity and interplanetary magnetic field vectors. In 
the direction across the magnetic field and the solar wind 
velocity, the shock results from the interaction of the plasma 
flow with the region of the enhanced magnetic field inside 
the cavity that plays the role of the magnetic barrier. 

[38] The other main features of the solar wind interaction 
with the Moon known from observations and other numeri- 
cal experiments are reproduced in this idealized ID hybrid 
model: the formation of the cavity in the solar wind plasma 
and its collapse with ion-acoustic velocity, propagation of 
the rarefaction wave, negative charging of the cavity, exis- 
tence of counterstreaming proton beams parallel to the 
magnetic field inside the cavity and in the rarefaction region, 
thus, validating our approach. 

[ 39 ] For solar wind parameters (7] cT e ~ 100 eV, i.e., hot 
solar wind streams) used in the simulation, the shock starts 
to form at the approximate distance ~3.5 R M . The distance 


grows as the ratio of the IMF components B v /B x decreases 
reaching 7 R M for B v /B x ~ 0.5 (further decrease of the ratio 
causes very fast growth of the cavity length). The size of the 
region occupied by the shock is rather small. The distance 
between the quasi-parallel shock and the wake axis is only 
0.4 R m atX= 9 R m . We repeated the simulations with lower 
electron temperatures (T e ~ 20 eV) and found weakened 
shock fonnation trailing the moon at much greater distances. 
Therefore, in order to observe the trailing shock, a satellite 
should have a trajectory passing very close to the wake axis 
during the period of hot solar wind streams. In principle, 
ARTEMIS spacecrafts placed at LI and L2 points of the 
Earth-Moon system [Angelopoulos, 2011] may cross the 
lunar wake sufficiently far downstream. 

[40] Acknowledgments. LO would like to acknowledge support by 
NASA grant NNX10AC56G. 
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